The automatic hp adaptation algorithm was first described in detail in the books of prof. Leszek Demkowicz, a Polish mathematician working at the University of Texas in Austin. The analysis of the mathematical properties of the hp-adaptation Algorithm was performed by prof. Ivo Babuška, an American of Czech descent, also working at the University of Texas in Austin [1], [2], [3], [4].

The hp-adaptation algorithm is illustrated in Fig. 1, Fig. 2, Fig. 3, Fig. 4, Fig. 5.
The convention for drawing computational meshes is as follows. Colors in the drawings Fig. 1, Fig. 3 and Fig. 5 denote the degrees of polynomials spanning over the edges and interior of the elements. One polynomial degree is provided on each edge, while a horizontal degree and a vertical degree are provided on each interior of the element.
For example, in the picture Fig. 1 all polynomials on all edges have degree 2, and all polynomials on the inside of all elements have degree 2 in the vertical and horizontal directions. In the picture Fig. 3 all polynomials on all edges have degree 3, and all polynomials on the inside of all elements have degree 3 in the vertical and horizontal directions. See the drawing Fig. 5. When viewing elements by rows from left to right, in the first row on the first element all polynomials on the edges have degree 2, except for the polynomial on the lower edge of the element which has degree 1. On the first element the polynomials in the interior have degree 2 in the vertical and horizontal directions. On the second element, on the left, top and bottom edges, the polynomials have a degree 2, while on the inside, a degree 3 in the vertical direction and a degree 2 in the horizontal direction. On the third and fourth elements, all polynomials on the edges are of degree 3. Likewise, the polynomials in the interior in the vertical and horizontal directions. In the second row of items, the first item (i.e., the fifth global) has degree 3 polynomials on the left and right edges, and a degree 1 polynomial on the top and bottom edges. Inside, the element has degree 1 polynomials in the horizontal direction and degree 3 polynomials in the vertical direction. The second element in the second line has a polynomial degree 2 in the horizontal direction and a polynomial degree 3 in the vertical direction. On the third and fourth elements, all polynomials on all edges have degree 3, and all polynomials on the inside of all elements have degree 3 in the vertical and horizontal directions. The remaining elements in the third and fourth lines have polynomials symmetrically distributed with the polynomials in the first and second lines on the first and second elements.
What is the optimal method of adapting a single element? Note that in the case of the hp adaptation algorithm, we have many possibilities to modify a single element:

1. Leaving the item unchanged
2. It is possible to break an element into two new elements vertically
3. It is possible to break an element into two new elements in the horizontal direction
4. It is possible to break an element into four new elements (two in the horizontal direction and two in the vertical direction).

Moreover, for each broken element it is possible to modify the degree of polynomials inside the element

1. Leave element orders unchanged
2. It is possible to raise the order inside the element horizontally \( (p_h,p_v)\rightarrow (p_h+1,p_v) \)
3. It is possible to raise the order inside the element in a vertical direction \( (p_h,p_v)\rightarrow (p_h,p_v+1) \)
4. It is possible to raise the order inside the element horizontally and vertically \( (p_h,p_v)\rightarrow (p_h+1,p_v+1) \)

The orders at the edges of the modified elements are determined on the basis of the minimum rule. In other words, the degree of the polynomial at the edge divided by two adjacent elements is set as a minimum of degrees inside the elements in the appropriate direction.
For a vertical edge surrounded by two elements \( (p^1_h,p^2_v) \) and \( (p^2_h,p^2_v) \) we set the degree \( p=\textrm{min} \{ p^1_v,p^2,v\} \).
For a horizontal edge surrounded by two elements \( (p^1_h,p^2_v) \) and \( (p^2_h,p^2_v) \) we set the degree \( p=\textrm{min} \{ p^1_h,p^2,h\} \).
So we have a lot of possibilities to adapt a single element (estimating on the basis of the possible types of adaptation we have 4 (modifications of the degree of the element without breaking the element) + 4 * 4 (breaking the element into two elements in the horizontal direction and four possibilities of modifying the degree of each element) + 4 * 4 (break an element into two elements vertically and four possibilities to modify the degree of each element) + 4 * 4 * 4 * 4 (break an element into four elements and four possibilities to modify the degree of each element). A total of 4 + 4 * 4 + 4 * 4 + 4 * 4 * 4 * 4 = 292 possibilities.
How is the decision about the type of adaptation of a single element made?
The decision is made in accordance with the following Algorithm.

In other words, we choose the type of adaptation for the element that gives the greatest error reduction at the lowest cost. This quantity is represented by the error decrease rate, which increases with the decrease in error but decreases with the cost incurred (with the number of functions added to the element because the cost of the calculation on a given element depends on the number of unknowns which are coefficients of the basis functions). This algorithm is illustrated in Fig. 6.

For the exemplary problem of heat transport in the L-shaped area, it is possible to obtain an adaptive hp mesh that allows solving the problem with the relative error of 0.0001. Such a grid is shown in Fig. 7, Fig. 8, Fig. 9, Fig. 10, Fig. 11, Fig. 12, Fig. 13.

The automatic hp adaptation algorithm is therefore very complicated and difficult to implement. So why is it worth our interest? Consider possible different adaptation algorithms for the exemplary L-shaped heat transport problem, i.e.

1. A homogeneous adaptation algorithm where the number of elements is increased evenly over the entire area
2. A homogeneous p adaptation algorithm, using a mesh made of 3 elements and increasing uniformly the degree of polynomials inside all elements horizontally and vertically, and the degree of polynomials on all edges
3. A homogeneous p adaptation algorithm, using a grid made of 12 elements and increasing uniformly the degree of polynomials inside all elements horizontally and vertically, and the degree of polynomials on all edges
4. An automatic p adaptation algorithm, increasing the degree of polynomials in the interiors of selected elements in selected directions, and modifying the degrees on the edges according to the minimum rule
5. An automatic adaptation algorithm, breaking selected elements in selected directions
6. An automatic hp adaptation algorithm, described in this chapter

If we draw a graph of the convergence of these algorithms in such a way that on the horizontal axis we place the size of the grid understood as the number of basis functions on the entire grid, which corresponds to the number of unknown coefficients of these basis functions, i.e. the size of the matrix of the system of equations to be solved, and on the vertical axis, the relative error calculated in H1 standard between the solution on the coarse mesh and the dense mesh
\( \|u_{h,p}-u_{h/2,p+1}\|_{H^1}=\frac{\int_{\Omega} (u_{h,p}-u_{h/2,p+1})^2+(\frac{\partial (u_{h,p}-u_{h/2,p+1}}{\partial x})^2+(\frac{\partial (u_{h,p}-u_{h/2,p+1}}{\partial x})^2 }{\int_{\Omega} (u_{h/2,p+1})^2+(\frac{\partial (u_{h/2,p+1}}{\partial x})^2+(\frac{\partial (u_{h/2,p+1}}{\partial x})^2 } \)
then it turns out that only the automatic hp adaptation algorithm has the feature of exponential convergence. It is definitely the fastest and can solve a given computational problem with an accuracy that is impossible to achieve by other adaptive algorithms.

then it turns out that only the automatic hp adaptation algorithm has the feature of exponential convergence. It is definitely the fastest and can solve a given computational problem with an accuracy that is impossible to achieve by other adaptive algorithms.
In practice, most engineering problems require accuracy in the range of 5 percent and square polynomials. However, there are computational tasks where high accuracy of the solution is essential. Examples of such computational tasks are, for example, simulations of flows around the wings of aircraft, which require very high accuracy in the boundary layer on the wings of the aircraft, or simulations of the propagation of electromagnetic waves in the geological formation layers to identify oil deposits that require very high accuracy in the vicinity of the receiver antenna.